function Ygroups=ssa1(Y,L,grouping,doplot)
%Computes groups of principal components of the time series Y corresponding to
%indices specified in the vector "grouping" which identifies the first index
%of each grouping of principal components. The lowest index corresponds to the
%highest variance (eigenvalue)


N=length(Y);
Noutputseries=length(grouping);

Ymean=mean(Y); Y=Y-Ymean;
X=zeros(L,N-L+1);
for k=1:N-L+1
    X(:,k)=Y(k:k+L-1);
end

[U,D] = eig(X*X'); D=diag(D); U=fliplr(U); D=flipud(D);

Ygroups=zeros(Noutputseries,N);

for n=1:Noutputseries-1
    if n<Noutputseries
        lastindex=grouping(n+1)-1;
    else lastindex=L;
    end
    
    %sum of elemantary matrices from projections onto eigenspace
    Xnew=zeros(size(X));
    for k=grouping(n):lastindex
        %[lastindex n k]
        Xnew=Xnew+U(:,k)*U(:,k)'*X;
    end

    %diagonal averaging
    Ynew=zeros(1,N);
    for j=1:N-L+1
        a=1; b=j; c=0;
        while a<=L & b>=1
            Ynew(j)=Ynew(j)+Xnew(a,b);
            a=a+1; b=b-1; c=c+1;
        end
        Ynew(j)=Ynew(j)/c;
    end
    for j=N-L+2:N
        a=j-(N-L); b=N-L+1; c=0;
        while a<=L & b>=1
            Ynew(j)=Ynew(j)+Xnew(a,b);
            a=a+1; b=b-1; c=c+1;
        end
        Ynew(j)=Ynew(j)/c;
    end

    %storing result
    Ygroups(n,:)=Ynew;
end
Ygroups(Noutputseries,:)=Y-sum(Ygroups(1:Noutputseries,:));
    
if doplot
    figure(1); clf;
    subplot(Noutputseries+1,1,1); plot(Y,'r'); ax=axis(gca);
    for n=1:Noutputseries
        subplot(Noutputseries+1,1,n+1); plot(Ygroups(n,:)); %axis(ax);
    end
    figure(2); clf; hold on; plot(D); plot(D,'r.'); set(gca,'yscal','log');
end